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Bose-Einstein Condensation Dynamics from the Numerical 
Solution of the Gross-Pitaevskii Equation by the 
Pseudospectral Method 



Paulsamy Muruganandam and Sadhan K. Adhikari 



Abstract — We study certain stationary and time-evolution problems of 
trapped Bose-Einstein condensates of weakly interacting alkali atoms de- 
scribed by a nonlinear Gross-Pitaevskii (GP) equation. We suggest a pseu- 
dospectral method involving Laguerre polynomials to solve the time depen- 
dent GP for a spherically symmetric trap potential. The radial wavefunc- 
tion and energy values have been calculated for different nonlinearities. 
Further, we study the effect of suddenly changing the interatomic scatter- 
ing length or harmonic oscillator trap potential in the condensate. We also 
investigate the frequency of oscillation due to the variation in the strength 
of nonlinearity. 

Keywords — Bose-Einstein condensation, Pseudospectral method 



I. Introduction 

THE experimental realization! 1 1 of Bose-Einstein conden- 
sates (BECs) in dilute weakly-interacting trapped bosonic 
atoms at ultra-low temperature initiated intense theoretical ef- 
fort to describe the properties of the condensate. The proper- 
ties of a condensate at zero temperature are usually described 
by the time-dependent, nonlinear, mean-field Gross-Pitaevskii 
(GP) equation| 1 1. The effect of the interatomic interaction (few- 
body correlation) leads to a nonlinear term in the GP equa- 
tion which complicates the solution procedure. Although there 
have been previous studies of the solution of the GP equation 
for stationary or time-independent problems, virtually all time- 
dependent studies in realistic cases, e.g., in three space dimen- 
sions, have employed approximate approaches rather than exact 
numerical solution of the GP equation. 

A numerical study of the time-dependent GP equation is of in- 
terest, as this can provide solution to many stationary and time- 
evolution problems. The time-independent GP equation yields 
only the solution of stationary problems. As our principal in- 
terest is in time evolution problems, we shall only consider the 
time-dependent GP equation in this paper. 

Here we propose a pseudospectral time-iteration method for 
the solution of the three-dimensional GP equation with spher- 
ically symmetric trap. In the pseudospectral method the un- 
known wave function is expanded in terms interpolating polyno- 
mials |2|. When this expansion is substituted into the GP equa- 
tion, the (space) differential operators operate on a set of known 
polynomials and generate a differentiation matrix operating on 
the unknown coefficients. Consequently, the time-dependent 
partial-differential nonlinear GP equation in space and time vari- 
ables is reduced to a set of coupled ordinary differential equa- 
tions (ODEs) in time which is solved by a fourth-order adaptive 
step-size controlled Runge-Kutta method |3| using successive 
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time iteration. Reinhardt and Clark |4| have employed pseu- 
dospectral methods for the solution of the GP equation, where a 
variable step forth order Runge-Kutta time propagator was used, 
as in the present work. A pseudospectral Fourier-sine basis was 
used in |4| for finite traps, and a corresponding complex pseu- 
dospectral basis was used for systems with periodic boundary 
conditions. 

Such a pseudospectral method is used in |5| with Her- 
mite polynomials as the interpolant for the case of completely 
anisotropic trap potential. However, when the system has ra- 
dial or axial symmetry the three dimensional problem reduces 
to one or two dimensions. In such situations, it will be more ad- 
vantageous to employ the integration in reduced dimensions as 
it saves considerable amount of computer memory and time. In 
this paper we considered the GP equation for spherically sym- 
metric trap potential. 

In section|n|we describe briefly the three-dimensional, time- 
dependent GP equation with spherically symmetric trap. The 
pseudospectral method is described in section||ll] In section HVl 
we report the numerical results for the wave function and energy 
for different nonlinearities as well as an account of our study on 
the effect of sudden changes in the interatomic scattering length 
or harmonic oscillator trap potential in the condensate. An anal- 
ysis on the frequency of oscillation with respect to variation of 
the strength of nonlinearity is also reported. Finally, in section 
[y]we present our conclusions. 

II. Nonlinear Gross-Pitaevskii Equation 

At zero temperature, the time-dependent Bose-Einstein con- 
densate wave function ^(r; r) at position r and time r may be 
described by the following mean-field nonlinear GP equation 1 1 1 
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V(r)+gN \^(r;r)\ 2 -in^- 
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*(r;r) = 0. 



(1) 



Here m is the mass and Nq the number of atoms in the conden- 
sate, g = Airti 2 a/m the strength of interatomic interaction, with 
a the atomic scattering length. The normalization condition of 
the wave function is J dr|W(r; r)| 2 = 1. 

The three-dimensional trap potential is given by V(r) — 



\m{u: 2 x x 2 



9 -? 



-), where uj x = cjq, w y , and uj z are 



the angular frequencies in the x, y and z directions, respectively, 
and r = (x,y,z) is the radial vector. 

A. Spherically Symmetric Case 

In the spherically symmetric trap, i.e., u> x = ui y = uj z = ui 
the tran ootential is dven bv V(r) = hmuinr 2 . where un is 
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the angular frequency and f the radial distance. The wave func- 
tion can be written as ^(r; r) = ip(r, r). After a transforma- 
tion of variables to dimensionless quantities defined by r = 

y/2f/l, t = tuiq, I = \J (fr/mujQ) and cf>(r,t) = ip(r,t)/r = 
ijj(r, r)(47r/ 3 ) 1 / 2 , the GP equation in this case becomes 



<p(r, t) 



d 



<p(r,t) = 0, 



(2) 



where M = N^a/l. The normalization condition for the wave 
function is 



/•OO 

/ I ip(r,t) | 2 dr = 2V2. 
Jo 



(3) 



III. PSEUDOSPECTRAL RUNGE-KUTTA (PSRK) METHOD 

The main idea of this method is that the spatial differentiation 
operators in equation (|2j are replaced by appropriate differential 
matrices which are derived from spectral or pseudospectral col- 
location. Then the time integration is carried out by any ordinary 
differential equation (ODE) solver. 

In the pseudospectral approximation, an unknown function 
f(x) is expanded in terms of weighted interpolants of the form 
[2|,|6| 



f(x) ~p n (x) 



(4) 



Here {;Ej}™_ is a set of distinct interpolation nodes, fj = 
f(xj), a(x) is a weight function, the interpolating functions 
{£j(:e)}™ = q satisfy £j(xk) = $jk (the Kronecker delta) and in- 
volve orthogonal polynomials of degree (n), so that f(xk) = 
Pn{xk), k = 0,1, ... ,n. One could use orthogonal polyno- 
mial such as, Chebyshev, Hermite, Laguerre and Legendre as 
the interpolant. Even non-polynomial interpolants like Fourier 
(spectral) expansion of the function in terms of periodic cosine 
and sine functions can be used in the case of periodic boundary. 
The choice of interpolant depends on the domain and nature of 
boundary conditions. 

Differentiating Q m times, the m^ 1 derivative of f(x) at the 
nodes is given by, 



f (m) (x) 



^ dx" 



a(x) 
a(xj) 



(5) 



with k = 1, . . . ,n. Here the derivative operator is represented 
by £)( m ) which is the differentiation matrix with elements, 
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(m) 



dx r 



a{x) 
a(xj) 



(6) 



Thus, the numerical differentiation of the function f(x), /' m ' = 
DWf is performed as the matrix- vector product where / and 
j( m ) ^ the vectors of the function and its derivative values 
evaluated at the nodes. 

In the present work we consider Laguerre polynomials L n (x) 
as the interpolating functions with the weight function a(x) = 
e~ 2 bx p(x) . The reason behind the choice of Laguerre polyno- 
mials is that thev are well defined in the interval x £ fO. 00) and 



satisfy the boundary conditions of the wave function of the GP 
equation |7|. The choice x = r 2 gives exact eigenfunctions of 
the GP equation (0 when the nonlinearity term J\f = 0. Actu- 
ally, the associated Laguerre polynomials are more appropriate 
as they are the radial eigenfunctions of linearized GP equations. 
However, from a numerical point of view the use of both of the 
polynomials yield more or less the same result. Also it has not 
been known whether there is a significant advantage of improve- 
ment in the results by considering associated Laguerre polyno- 
mials^!. 

Consequently, the partial differential equation (0 is reduced 
to a set of coupled ODEs in the time variable t involving £j , j = 
0,1, ... ,n. In this way we obtain a set of ODEs by collocating 
the original equations on a suitable set of points (the roots of 
some Laguerre polynomials). 

For solving the set of ODEs we use the adaptive step-size 
control based on the embedded Runge-Kutta formulas due to 
Fehlberg |3|, which gives a definite clue about how to modify 
the step size in order to achieve a desired accuracy in a con- 
trolled way. For orders M higher than four of the Runge-Kutta 
formula, evaluation of more than M functions (though never 
more than M + 2) is required. This makes the classic fourth 
order method requiring the evaluation of four functions more 
economic. Fehlberg suggested a fourth-order and a fifth-order 
method each requiring the evaluation of six functions. The dif- 
ference between the results of these two gives the error S in the 
fourth-order method with a step size h, where S scales as h 5 : 
S oc h 5 . This scaling immediately gives the factor by which the 
step size h should be reduced, so that a desired S can be ob- 
tained. The detailed fourth-order and fifth-order Runge-Kutta 
formulas of Fehlberg are given in |3 |. We use these formulas 
with the constants given by Cash and Karp also tabulated in 1 3 1. 
For the present problem we find that the use of Cash-Karp con- 
stants in the Fehlberg formulas leads to more accurate results 
than the original constants due to Fehlberg. 

Using the differentiation matrix ©, the GP equation is dis- 
cretized. The grid points are the roots of the Laguerre polyno- 
mial L n (xj) = 0. However, the actual Xj values employed are 
obtained by scaling these roots by a constant factor so that most 
of the roots fall in the region where the condensate wave func- 
tion is sizable and only a few points are located in the region 
where the wave function is negligible. In the present paper, all 
the numerical simulations have been carried out with n — 31 
grid points. 

IV. Numerical simulation 
A. Wavefunction and energy 

Once the GP equation (0 is reduced to a set of coupled ODEs, 
the time integration is carried out with an initial seed wave- 
function. We use the solution corresponding to the linearized 
GP equation as the seed. The final wavefunction is obtained 
by adding the nonlinear terms in small steps until the desired 
nonlinearity is achieved. Then the stationary wave function <f> 
and the parametric energy /1 (the chemical potential) can be ex- 
tracted from the evolution of the time-dependent GP equation 
over a macroscopic interval of time Isl. l9l. 

First, we are interested to calculate the wavefunction and en- 
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ergy with different nonlinearity strengths. For this purpose, we 
numerically integrate equation with the following normal- 
ized analytic solution, 

?(r) = (4) r exp (-^) . (7) 

The above ip{r) corresponds to the ground state of (0 for spher- 
ical symmetry with the coefficient of the nonlinear term set to 
zero. During the integration the coefficient of the nonlinear 
term is increased from at each step by S n — 0.001 until the 
final value of nonlinearity J\f is attained at a time called time 
t = 0. This corresponds to the final solution and this solu- 
tion is found to be stable. The norm of the wavefunction is 
conserved during the integration. However, it is of advantage 
to reinforce numerically the proper normalization of the wave 
function given by (0 after finite number of RK steps in order 
to improve the precision of the result. As the grid points are 
the corresponding Laguerre roots, the integral in can easily 
be evaluated by appropriate Gauss-Laguerre formula. Figure 
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Fig. 1 

Figure showing the radial wavefunctions of the condensate 

FOR DIFFERENT NUMBER OF ATOMS Nq = 2«. 

depicts the radial wavefunction of the condensate for different 
nonlinearity strengths. The stationary wavefunction has trivial 



TABLE I 

THE CHEMICAL POTENTIAL fl FOR VARIOUS NUMBER OF CONDENSATE 
ATOMS. 



No 


M 


A* 


N 


M 


A* 


1024 


0.50 


1.824 


2048 


1.00 


2.065 


4096 


1.98 


2.429 


8192 


3.97 


2.969 


16384 


7.93 


3.722 


32768 


15.86 


5.748 


65536 


31.72 


6.155 


131072 


63.44 


8.041 



time dependence </>(r, t) = </>(r) exp(— ifxt). Thus, the paramet- 
ric energy /i (chemical potential) can be extracted from the evo- 
lution of the GP equation over a macroscopic interval of time. 
In Table [[] we present the chemical potential of several BECs 
with different number of atoms for the spherically symmetric 



trap potential. The typical value of the harmonic trap frequency 
is ujo = 87rad/s and scattering length is a = 52ao (for Na) 
where a is the Bohr radius|l|. The energy values calculated 
here for different nonlinearity strengths are compared reason- 
ably well with the results of Schneider and Feder 1 10 1 by time- 
independent approach. 

B. Time evolution 

Next we consider certain time-evolution problems which can 
be tackled more effectively with the present pseudospectral 
method. From the experimental point of view, it is of great inter- 
est to study the effects of changing (i) the trap potential and (ii) 
the atomic scattering length. The former can be achieved in the 
laboratory by changing the magnetic field responsible for trap- 
ping while it has been possible to modify the scattering length 
in a controlled fashion by exploiting a Feshbach resonance 11 II . 
In this paper we restrict ourselves to the study on the effect of 
suddenly changing both the trap potential as well as the atomic 
scattering length. 

B.l Effect of changing trap potential and scattering length 

As said above, one can increase or reduce suddenly the 
strength of the nonlinear term through a change of the scattering 
length in this fashion and study the oscillation of the condensate 
thereafter. In both these cases one observes the time evolution of 
the root mean square (rms) radius (r)rms of the condensate. For 
the time evolution study we consider a previously formed con- 
densate with finite number of atoms, for example, Nq — 16384 
prepared at t — as in the study of the stationary problem. We 
then inflict the four following changes in the system. At t = 
we (a) increase or (b) decrease suddenly the coefficient of the 
harmonic oscillator r 2 /4 term in (0 by a factor of 2. Next at 
t = we (c) increase or (d) decrease suddenly the coefficient of 
the nonlinear term JV in (|2} by a factor of 2. 

When the harmonic oscillator term is doubled or the nonlin- 
earity halved, the system is compressed and the rms radius os- 
cillates between its initial value and a smaller final value. When 
the harmonic oscillator term is halved or the nonlinearity dou- 
bled, the system expands and the rms radius oscillates between 
its initial value and a larger final value. 

We have calculated the frequency of oscillation of the rms 
radius in all the four cases considered above. These frequen- 
cies are calculated as, for the (a), (b), (c), and (d) cases above, 
0.4834, 0.2441, 0.3467, and 0.3418, (in units of lu^ 1 ) respec- 
tively. Here the time is measured in units of ui^ 1 = 1 and hence 
the trap frequency corresponding to Q is Vq = 1/(2tt). The 
frequency of oscillation was two times the existing harmonic 
oscillator frequency for nonlinearity Af = 0. As the frequencies 
are increased and reduced by \[2 in cases (a) and (b), two times 
the existing frequency for (a) is 2v = y/2/tr ~ 0.45 and for (b) 
is 2v = 1/(V2tt) ~ 0.23. In case of (c) and (d) the frequency 
is unchanged and 2v = 1/tt ~ 0.32. These numbers compare 
well with the respective above results. 

We also found that these frequencies are in good agreement 
with those calculated using finite difference scheme 1 12 1. There 
are differences when comparing these frequencies with the cor- 
responding harmonic trap frequencies, the reason for which is 
due to the presence of nonlinearitv in the svstem. We studv the 
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frequency of oscillation as a function of nonlinearity in the fol- 
lowing. 

B.2 Frequency of oscillation of the condensate 

Having studied various time independent as well as time de- 
pendent problems, we further investigate the frequency of os- 
cillation of the condensate due to a variation of nonlinearity. 
When the nonlinear term is added in sufficiently large steps dur- 
ing the preparation of the condensate numerically, a breathing 
type oscillation with some frequency, v occurs which is deter- 
mined solely by the final nonlinearity. This frequency can be 
extracted by exploiting the Fourier transform on the time series 
of the mean distance, (r) rms from numerical simulation. We 



0.36 




Fig. 2 

Plot showing frequency of oscillation versus the No 



express the final nonlinearity in terms of the final value of Nq. 
We found that the frequency of oscillation increases as the fi- 
nal number of atoms No is increased and it saturate to a final 
value 0.354 when N > 2 16 ~ 7 x 10 4 . Figure shows the of 
frequency oscillation v as a function of nonlinearity in the con- 
densate. By fitting the numerical data for the frequencies, the 
frequency of oscillation as a function of final Nq can be repre- 
sented by the following relation, 



u{N ) 



1 



■ exp 



1 



(8) 



-28.6858, 



The parameters are calculated as v' Q — 0.3536, b 
c = 227.26 and 7 = 0.6202. 

V. Summary and Conclusions 

We propose the use of pseudospectral method to solve time 
dependent nonlinear Gross-Pitaevskii equation describing Bose- 
Einstein condensation in spherically symmetric trap potential. 
In particular, we emphasize on the advantages of using Laguerre 
polynomial interpolation for the spatial differentiation operator 
and we study both stationary as well as time evolution prob- 
lems. In the stationary problem, the radial wavefunction and 
the ground state energy are calculated by solving the GP equa- 
tion for different nonlinearitv strengths which determines the 



number of atoms in the condensate. We noted that the energy 
values are calculated more accurately with very few number of 
discretization points (nodes). For the time evolution problem, 
first we study the effect of sudden changes in the interatomic 
scattering length or harmonic oscillator trap potential in the con- 
densate. We have also made an analysis on the variation of fre- 
quency of oscillation as a function of the strength of nonlinear- 
ity. 

This paper summarizes our attempts to understand certain dy- 
namical aspects of Bose-Einstein condensation as a nonlinear 
dynamical system. We found that the present approach based 
on pseudospectral differentiation matrices for solving Gross- 
Pitaevskii equation is more effective when compared to finite- 
difference schemes. Though, we are solving the problem effec- 
tively in one-dimension because of the spherically symmetric 
trap, the anisotropic problems can be easily tackled with similar 
approach(5). 
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